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ABSTRACT 

Gravitational instability of a vertically thin, dusty sheet near the midplane of a protoplanetary disk 
has long been proposed as a way of forming planetesimals. Before Roche densities can be achieved, 
however, the dust-rich layer, sandwiched from above and below by more slowly rotating dust-poor gas, 
threatens to overturn and mix by the Kelvin-Helniholtz instability (KHI). Whether such a threat is 
real has never been demonstrated: the Richardson criterion for the KHI is derived for 2-D Cartesian 
shear flow and does not account for rotational forces. Here we present 3-D numerical simulations of 
gas-dust mixtures in a shearing box, accounting for the full suite of disk-related forces: the Coriolis and 
centrifugal forces, and radial tidal gravity. Dust particles are assumed small enough to be perfectly 
entrained in gas; the two fluids share the same velocity field but obey separate continuity equations. We 
find that the Richardson number Ri does not alone determine stability. The critical value of Ri below 
which the dust layer overturns and mixes depends on the height-integrated metallicity Yid/Tig (surface 
density ratio of dust to gas). Nevertheless, for TjdfSg between one and five times solar, the critical 
Ri is nearly constant at ~0.1. Keplerian radial shear stabilizes those modes that would otherwise 
disrupt the layer at large Ri. If the height-integrated metallicity is at least ^5 times greater than the 
solar value of 0.01, then midplane dust densities can approach Roche densities. Such an environment 
might be expected to produce gas giant planets having similarly super-solar metallicities. 

Subject headings: hydrodynamics — instabilities — planets and satellites: formation — planetary 
systems: protoplanetary disks — turbulence 



1. INTRODUCTION 

How do dust grains, known to permeate disks surrounding young stars, assemble into planets? One proposed stage 
of growth involves gravitational instability of a dust-rich layer at the disk midplane (Safronov 1969; Goldreich & Ward 
1973). If dust grains are free to settle vertically out of gas, midplane dust densities eventually exceed Roche densities, 
and dust particles can begin to aggregate by self-gravity. 

An objection to this means of forming planetesimals is that vertical velocities of gas may be too high for dust to settle 
] (Weidenschilling 1980). Even apart from turbulence intrinsic to gas (e.g., turbulence driven by the magneto-rotational 
instability), dust layers that are too vertically thin can suffer from the Kelvin-Helmholtz instability (KHI). The KHI 
threatens to manifest because dust-rich gas at the midplane rotates at a rate different from that of dust-poor gas at 
higher altitude. Dust-poor gas experiences greater acceleration from a background radial pressure gradient, and so 
its rotation velocity in centrifugal balance must deviate more strongly from the Keplerian value. ^ The deviation is 
smaller for gas laden with dust, since dust adds inertia but contributes no pressure. Usually it is assumed that the 
background radial pressure gradient dP/dr < so that the vertical shear in the rotation velocity dv^/dz < 0. 

Whether the disk is KH- unstable is commonly assessed using the Richardson number (Chandrasekhar 1961, page 
■ 491): 

^. _ {g/p)dp/dz 
~ {dv^/dzY 

where g is the vertical gravitational acceleration and p is the total mass density. In the limit that perturbations 
are incompressible, the numerator is the square of the Brunt- Vaisala frequency of buoyant oscillations, while the 
denominator is the square of the vertical shearing rate. The Richardson number measures the amount of work 
required to overturn fluid elements that are originally in hydrostatic equilibrium, against the amount of free kinetic 
energy available in the background shear. For purely Cartesian flow (no rotational forces), 

Ri < Riciit = 1/4 is necessary for instability (1) 

(Miles 1961; Howard 1961; Chandrasekhar 1961; Drazin & Reid 1981; Li et al. 2003). Despite the fact that (1) is 
formally not a sufficient criterion, laboratory experiments bear out its usefulness (Tritton 1988). 
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In the case of our disk, close to the midplane, g = — fl^z arises from the tidal gravity of the star, z measures distance 
from the midplane, is the Keplcrian angular frequency, and p = Pg + Pd, where g and d denote gas and dust, 
respectively. Moreover, \dpg/dz\ <C \dpd/dz\ within the thin dust layers of interest to us (see §2.2). Thus throughout 
this paper 

Pd + Pg {dv^/dzf ■ 

As dust settles vertically, pd ~ T,d/ Zd increases; the dust surface density is assumed constant while the characteristic 
scale height Zd of dust decreases. In the limit that pd ^ Pg, Ri oc z^. Criterion (1) for instability therefore places a 
lower bound on Zd and a corresponding upper bound on pd- Unfortunately, for conventional solar nebular parameters — 
in particular for a height-integrated solar metallicity of Tid/'^g = 0.01— the maximum of pd falls short of the Roche 
density by 1 2 orders of magnitude (e.g., Sekiya 1998; Youdin & Shu 2002). 

But is it appropriate to apply criterion (1) to a circumstellar disk where the flow is not Cartesian? One might 
guess that rotational forces merely modify Riait by a factor of order unity, since the Coriolis force and the radial 
shear {(Kl/dr) operate on a timescale — the same timescale characterizing vertical Brunt- Vaisala oscillations and 
the vertical shear at z ^ Zd- Yet when Coriolis forces are introduced, as in the analysis of Gomez & Ostriker (2005, 
hereafter GO), the stability properties of the flow change dramatically. Numerical simulations by GO, performed 
under the assumption that dust is perfectly entrained in gas, reveal no well-defined threshold for instability. Unstable 
modes are detected for Ri as high as 4, with growth rates that diminish with increasing Ri but that show no sign of 
vanishing (see their Figures 9 and 10). Similar results obtain in numerical simulations by Johansen, Henning, & Klahr 
(2006, hereafter JHK), who relax the assumption of perfect coupling between dust and gas but who, like GO, retain 
only Coriolis forces and ignore radial shear. For their "rocks" and "pebbles" with momentum stopping times 10-50 x 
shorter than il^^, the dust distribution evolves to one where Ri increases from ^1 at the midplane to > 10 at higher 
altitude (JHK, their Figures 4 and 5). 

Here we conduct numerical simulations that account for the full complement of disk-related forces. We include not 
only the Coriolis force but also the centrifugal force and radial tidal gravity. The latter two forces combine to produce 
radial shear in rotational equilibrium. Thus we investigate the stability of doubly shearing flows: vertically shearing 
dust layers in near-Kcplcrian differential rotation. We expect our results to differ from those of GO and ,JHK. The 
modes emphasized by GO have growth rates < O.lfi at Ri > 2. Keplerian differential rotation, characterized by a 
strain rate of 30/2, should shear such modes apart before they have time to amplify. 

Indeed this is the conclusion of Ishitsu & Sekiya (2003, hereafter IS), who numerically integrate the linearized 
equations of motion, including the full suite of disk-related forces, to compute the factors by which modes amplify. 
They find that amplification factors are sufficiently limited by radial shear that midplane dust-to-gas density ratios 
might reach values as high as ~2 (still too low for gravitational instability, unfortunately). They restrict, however, 
their linear analysis to odd-parity modes for which the vertical velocity Vz{z) = —Vz{—z). Gomez & Ostriker (2005) 
establish that even-parity modes grow faster and overturn the dust layer more effectively. Our simulations do not pre- 
select for either type of mode. Another difference between our work and IS is that we concentrate on Ri = constant 
flows, whereas IS employ a simple, analytically tractable dust distribution. Of course, Ri = constant is merely a 
plausible condition to which the dust-gas mixture might relax. One of the goals of this study is to assess whether the 
Richardson number is a good predictor of stability even when all disk-related forces are accounted for. Finally, the 
analysis of IS is linear, while our numerical simulations enable access to nonlinear phenomena. 

Our numerical simulations, like those of GO, assume that dust particles are perfectly entrained in gas, i.e., we assume 
that grain momentum stopping times <C fi^^. This approximation is valid for small particles, e.g., having sizes <C 1 m 
at disk radius r = 1 AU in a minimum- mass nebula (e.g., Weidenschilling 1977). The hydrodynamics code and initial 
conditions are described in §2. Results are presented in §3; there we determine Ricvit for our doubly shearing flows. A 
summary is provided in §4. 

While our entire study is rooted in the literature on the Kelvin-Helmholtz instability and its impact on planetesimal 
formation, and as such revolves around the Richardson criterion, it was brought to our attention by the referee that the 
KHI may not be the only instability afflicting the dusty midplane. Baroclinic instabilities, studied in disks by Cabot 
(1984), Knobloch & Spruit (1986), and Arlt & Urpin (2004) may also be relevant. These authors study baroclinity 
specifically in the context of a non-zero vertical shear. More generally, a baroclinic flow is one whose isobaric surfaces 
do not coincide with its isodensity surfaces. Our flows are baroclinic because of molecular weight gradients: dust adds 
to the density but not to the pressure. Baroclinic instability can strike even when the Richardson number is large. We 
will connect our findings to the baroclinic instability when appropriate; indeed, in §3.1 we propose that the instability 
discovered by GO in their simulations of non-radially-shearing disks is, in fact, the baroclinic instability discovered 
by Cabot (1984). And we will see in §3.2.2 that the Richardson criterion alone does not determine stability under 
arbitrary conditions. 

2. METHOD 

2.1. Dust-Gas Equations in Tight Coupling Limit 
Gas and pressureless dust are two fluids obeying separate momentum equations. In an inertial frame, 
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^ + K^vK.-v* + ^i&pil-^ (2) 

CI Pg Cstop Pg 

^ + (v..V)v. = -V^-(^ (3) 

where v is velocity, t is time, and $ is the background potential. Terms containing tstop, the momentum stopping time 
for dust in gas, account for how the two fluids interact frictionally. 

In this paper, we work exclusively in the tight coupling limit where particles and gas have negligible relative velocity, 
|vd — Vgl ^ |vg|. This corresponds to the case where dust particles are so small — i.e., their surface area to mass 
ratios are so large — that they become entrained in the gas flow over short times. By computing the difference and 
appropriately weighted sum of (2) and (3) (see Youdin & Goodman 2005), and taking the limit <stop 0, we find 

vd-vg = (4) 
^ + (vV)v = -V$-^— (5) 

5* Pp + Pg 

where v = = Vg. 
Furthermore, 

^ = - r^—^ (6) 
P=(7-l)e (7) 

de 

_+V.(£v) = -PV-v (8) 

where G is the gravitational constant, M is the central stellar mass at the origin, e is the internal energy density of 
gas, and r and z are the cylindrical radius and height. The gas obeys a polytropic equation of state P = Kp^g with 
constant K and 7. Finally, dust and gas obey separate continuity relations 

^ + V-(pgv)=0 (9) 

^+V.(p,v)=0. (10) 

Equations (5)-(10) are the ones we solve in this paper, in the shearing box approximation (Goldreich & Lynden-Bell 
1965; Hawley, Gammie, & Balbus 1995; §2.2). They contain the same content as equations (2)-(6) of GO. 

In the tightly coupled limit, a fluid parcel does not change its dust-to-gas ratio in a Lagrangian sense: d{pd/ Pg)/dt = 
0, where d/dt is the convc!ctive derivative. Particles cannot slip away from gas, and so we cannot expect our simulated 
flows to relax, from arbitrary initial conditions, into a unique, steady state of marginal stability (if such a state actually 
exists) . In reality, such relaxation occurs over the time it takes dust to settle vertically relative to gas. Vertical settling 
times are measured in Myrs for micron-sized particles (or ~10^ yr for centimeter-sized particles) in a minimum-mass 
nebula. Nevertheless, our equations do permit dust-rich parcels to settle toward the midplane if their weight cannot be 
supported in vertical hydrostatic equilibrium. We are able to evolve dust-gas mixtures over dynamical times, measured 
in orbital periods, provided tgtop is still shorter. Our simulations can thus determine whether a given set of equilibrium 
conditions is dynamically stable. We now specify these equilibrium initial conditions. 

2.2. Initial Conditions: Constant Ri Flows 



Call ro and (t>Q — ^.^t the radius and azimuth of a test particle moving on a circular orbit, where fio = \/ GM/tq. 
Shifting to axes centered on the test particle and rotating at Oqi we trade the usual cylindrical coordinates (r, </>, z) 
for their shearing sheet counterparts {x,y,z): x = r — tq, \x\ <C tq, y = (^ — (t>Q)rQ, \y\ 4C ro, and \z\ <C ro- In this 
rotating frame, the momentum equation (5) reads 



dvx dvx —1 dP 
dt * dxi Pg + pd dx 
dvy dv„ —1 dP 



+ 2noVy + 2qn^x (11) 



+ Vi^ = : ^-2QoVx (12) 



dt dxi Pg + Pd dy 
dvz dvz —1 dP 
dt ' dxi Pg + Pd dz 



nlz (13) 



4 



Chiang 



whore i = (1.2,3), (xi, X2, x^^) = {x, y,z), (wi,?^2,i'3) — {vx,Vy,Vz), and curvature terms of order u^/r are dropped, 
following the shearing sheet approximation. We introduce the variable q = —dlnQ/dhir in (11) because in our 
simulations we experiment with q = (zero radial shear) to connect to GO. For the full problem with a point-mass 
potential, q = 3/2. 

In rotational equilibrium, equation (11) for radial momentum balance gives 

1 dP 

Vy = -qflox + TTTT— ■ r -jr- ■ (14) 

2no(pd + Pg) OX 

Initially, dP/dx = {dP/dr)t=o reflects how the background nebular pressure changes radially on scales of r. We 
express {dP/dr)t=Q in terms of a model input parameter, Wmax: 



dP 
dr 



= -2pg{z, t = 0)i;maxfio • (15) 

t=0 



The velocity Wmax represents the difference between the Kcplcrian rotation rate and the sub-Keplerian rotation rate 
of pressure-supported gas; it measures the maximum possible difference in rotational velocity between the dust-rich 
midplane and dust-free gas at higher altitude. For typical nebular parameters, it is approximately 20 m/s and nearly 
constant with r. 

Further defining the local dust-to-gas ratio /x = Pd/Pg, we rewrite (14) as 

Vy{x,z,t = 0) = -qQoX- Y^^^- (16) 

In each of our simulations, we hold f,nax fixed for simplicity. Then the vertical shear, dvy/dz, is non-zero only when 
dji/dz is non-zero. Because /i varies rapidly with z for the vertically thin dust layers of interest to us, approximating 
^^max as constant introduces negligible error. 
All our simulations investigate the stability of constant Ri flows. The condition Ri = constant yields /u(z): 



Ri = -- 



Ctlz dpd/dz 

i + Pg [dvy/dzf 



»^ (1 + M)^^ 
<ax c)^^l^z 

which integrates to 



li{z) = 



1 



ll{l + ^QY + {zlzd) 
where iiq = ^{z — Q) is a model input parameter and 

Ri^/'^v 

is a characteristic dust height. Equation (17) implies that dust extends to a maximum height 

_ /xy^(2 + /xo)V^ 

2max — Zd ^ 

1 + /U0 

Finally, vertical hydrostatic equilibrium gives Pg{z): 

1 dP 



1/2 

- 1 (17) 



Pg + Pd dz 



n^z (18) 



which, using (17), can be solved analytically for pg. The expressions, which differ for < 2; < Zmax and z > ^max, are 
not especially illuminating and so we omit them here. 

Self-gravity is ignored. Its neglect is justified for densities less than the Roche density, {pg + pd) <C Af/(27rr'^). 
Equivalently, for standard nebular parameters, we are restricted to p <^ 30. Sekiya (1998; see also Youdin & Shu 2002) 
relax;es this condition and accounts for vertical self-gravity in computing the Richardson number. 

Figure 1 shows profiles for Pg{z) and Pd{z), for normalized (code) parameters of pgQ = Pg{z = 0) = 1, if = 1, 

7 = 5/3, f^o = 1, -Ri = 1, Pa = 0.475, and fmax/cso — 0.05, where c^o = \J iKp^^^ = 1.3 is the sound speed at 

the midplane. These parameter choices imply a height-integrated dust-to-gas surface density ratio (height-integrated 
metallicity) of S^/Eg = 0.01, the nominal solar value. We vary Ri and po from simulation to simulation. Note that 
our choice of Vmax/cso = 0.05 is a factor of 2 smaller than that of GO because we prefer to use gas temperatures 
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Fig. 1. — Dust and gas vertical density profiles, computed for constant Ri = 1 and vertical hydrostatic equilibrium. Disk self-gravity is 
neglected (cf. Sekiya 1998). Note the change in horizontal scale between the two panels. The dust density is zero at z = Zmax. Typically 
our simulations span z = ±2zmiix', the gas density pg is practically constant on these scales. 

characterizing passive disks (see, e.g., Chiang et al. 2001) rather than the hotter temperatures of the Hayashi (1981) 
nebula. Disk gas must be fairly passive, i.e., not heated by turbulent dissipation, for dust to settle into the thin layers 

of interest here (see Appendix B of Youdin & Chiang 2004 for estimates of the degree of turbulence permitted) . A 
factor of 2 decrease in Wmax/cso produces a factor of 2 increase in the central concentration /io, at fixed Ri and fixed 
T,^/Tig. Figure 2 is analogous to Figure 1, except that Ri = 0.125 and /ig = 1-31 (so that again Sj;/Eg = 0.01). 
Figure 3 is analogous to Figure 2, except that /xq = 23.3, so that Sd/^s = 0.05 (the height-integrated metallicity is 
super-solar by a factor of 5). 

To summarize this subsection, our equilibrium initial conditions specify Vy through (16), fi = Pd/Pg through (17), 
and Pg through analytic solution of (18). Initially, Vx = Vz = 0, and pg and pd are constant with x and y. The primary 
input parameters are Ri and /io; we hold Wmax/cso = 0.05 fixed for all runs. 

2.3. Code 

We employ the hydrodynamics code ZEUS (Stone & Norman 1992), a version of which was kindly given to us by 
Dr. Bryan Johnson. Because we are interested in flows that shear both in radius {dvy/dx) and height (dvy/dz), the 
simulations are necessarily three-dimensional. We adopt the usual shearing box (e.g., Hawley, Gammie, & Balbus 
1995), with shearing periodic boundary conditions in x, periodic boundary conditions in y, and closed boundary 

conditions in z. 

Three modifications are made to the code. The first is the addition of another fluid, dust, which is transported 
according to its own continuity equation (10), but otherwise has the same velocity as that of gas. The second 
modification is to the pressure gradient source term: VP/pg — > VP/{pg + pd), in accordance with (2.2). The final 
change is to introduce an extra source term for the "large-scale" radial pressure gradient, distinct from the just 
mentioned "local" pressure gradient. That is, equation (11) for the x-momentum is revised to read 

dvx dvx -1 dP „^ „ ^2 1 fdP\ 

-j^+Vi-^ = ^—— + 2^QVy + 2qnlx ■ 19 

dt dXi pg +Pd ox Pg+Pd J 

The extra source term {dP/dr)t=o is given by (15). It accounts for how the pressure changes over radial lengthscales 
of r, and is held fixed for each simulation. Its inclusion is necessary to have dust-rich gas rotate faster than dust-poor 
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z (normalized code units) 



Fig. 2. — Same as Figure 1, but for Ri = 0.125. 

gas initially, i.e., to drive the vertical shear. But because our simulations use a shearing box in other words, because 
they are local — we cannot follow the evolution of this large-scale gradient. A self-consistent treatment would require 
us to model the disk globally. Our use of an external, constant, radial pressure gradient is the same device employed 
by GO. The error introduced is negligible insofar as changes in pressure over lengthscales of r occur over timescales 
much longer than the dynamical timescales which concern us here. 

In short, the modified code solves equations (7) (10), (12), (13), and (19). Initial conditions for the code are given 
above in §2.2. The equilibrium state is perturbed by adding a random velocity to each grid cell. The perturbation 
velocity points in a random direction from cell to cell. For our main set of runs, the magnitude of the perturbation 
in each cell is drawn at random from a uniform distribution that extends from zero to 10~^Cso (for comparison, the 
maximum vertical shearing velocity is about Wmax = O.OScso)- Several auxiliary runs use much smaller or more spatially 
limited initial perturbations; these arc described in §3.2.1 and produce the same outcomes as our main simulations. 

Our standard box dimensions are = Ly = S^max and = 4:Znia.x, and the corresponding number of grid cells is 
{Nx, Ny, Nz) = (64, 64, 32). The duration of each simulation is at least tf — 20^^^. Run parameters are contained in 
Tables 1 and 2, and are justified in §3. 

Including the von Neumann & Richtmyer (1950) artificial viscosity makes little difference to our results — not sur- 
prisingly since our flows are highly subsonic (t^max/cso = 0.05). Therefore the simulations reported in this paper have 
no artificially imposed viscosity apart from the unavoidable numerical kind. 

3. RESULTS 
3.1. Coriolis Only: q = 

To test our code, we repeat the experiments of GO. The Coriolis acceleration is retained but the radial shear is 
suppressed by setting q = 10~^ in the code. Effectively, this converts shearing boundary conditions in x to periodic 
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Fig. 3. — Same as Figure 2, but for a super-solar metallicity S<j/Sg = 0.05. The Richardson number, 0.125, is about the lowest it can 
be without the dust layer turning over for this metallicity (see Run S8 of Table 2). The midplane dust density, more than 20 times higher 
than the local gas density, is sufficiently high that dust self-gravity can be dynamically important. Ways of achieving S^/Sg = 0.05 are 
mentioned in §4. 

boundary conditions. It also eliminates the contribution, 2gfi§a;, to dvx/dt from the centrifugal force and tidal gravity. 
Since the resultant flows have no structure in x, we reduce iVx to 2 and Lx to 0.25^;max- 

Table 1 lists the run parameters. The main parameter varied from run to run is Ri; the central concentration /io 
is adjusted in tandem to keep Ti^fSg = 0.01 fixed for all runs (this restriction to solar metallicity is relaxed for the 
fully shearing runs of §3.2). Following GO, we diagnose the flow by Fourier analyzing in the y-direction at fixed 
X, z, and t. We then inspect how the Fourier amplitudes change with time. We are interested in those modes whose 
y- wavelengths Xy are comparable to the dust layer thickness Zmax, as shorter wavelength modes cannot overturn the 
layer and longer wavelength modes grow more slowly. We are able to measure the exponential growth rates toi of 
several modes whose Aj,'s range from ^8/15 to ^8/3 of Zmax.^ Figure 4 displays the Fourier amplitudes versus time 
for the Xy/Zmax = 8/7 mode, in simulations of varying Ri. While the growth rate loj of this mode decreases with 
increasing Ri, growth is still reliably detected for Ri as high as 16, a value 64x greater than the canonical threshold 
of 1/4. Our measured growth rates are similar to those reported by GO (see their Figure 10). 

Runs Cla-Clf test the sensitivity of our results to simulation box size and spatial resolution. Doubling the box 
height Lz to Sz^ax at fixed resolution produces negligible change in measured growth rates (Cla vs. CI). Doubling 
the vertical resolution by doubling Nz at fixed Lz increases mode growth rates by ^^10% (Clb vs. CI). Changes in 
u)i of order 10% are also produced by varying the azimuthal parameters Ly and Ny by factors of 2 (Clc Clf). We 
conclude that our standard box size {Ly,Lz) = {8,4:)zmax and resolution {Ny,Nz) = (64,32) adequately balance the 
need for accuracy with the need for computational speed. 

Table 1 and Figure 4 suggest that when only the Coriolis force is included, and radial shear is omitted, dust layers 
characterized by surprisingly large values for Ri will eventually overturn and mix. This conclusion is supported by 
Figure 5, which displays snapshots of a high-resolution {Ny,Nz) = (128,64) run for which Ri = 4. Despite Ri being 
16 X greater than the traditional critical value of 1/4, the dust layer is clearly unstable. The fastest growing mode 
has Ay/zmax ~ 8/5 and overturns the layer by f w lOO^o ^. All these results agree with those of GO. 

We suspect that the instability discussed in this section, and simulated first by GO, is a baroclinic instability aflilicting 

^ Periodic boundary conditions in y imply that the only modes present are those for which the simulation box fits an integral number of 
wavelengths. The standard box length is Ly = Szmax. 
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TABLE 1 
CoRiOLis Only {q = 0) Simulations 



Name 




{•2max) 


Ny 


Ly 

(^^max) 




(■S^max) 


Ri 


Mo 




*^ 


Ay 

(•2max) 


(Ho) 




2 


0.25 


64 


8 


64 


8 


0.25 


0.903 


0.01 


20 


8/3 


0.334 


Clb 


55 


)5 


64 


8 


64 


4 




55 


55 


55 


8/3 


0.371 


Clc 


55 


)5 


64 


16/3 


32 


4 




55 


55 


55 


8/3 


0.352 


Old 


55 


)5 


32 


8/3 


)5 


55 




55 


55 


55 


8/3 


0.368 


Cle 


55 


)5 


32 


3 


55 


55 




55 


55 


55 


3 


0.339 


Clf 


55 


55 


64 


7 


55 


55 




55 


55 


55 


7/3 


0.336 


01 


55 


)5 


64 


8 


32 


4 




55 


55 


40 


8/3 
8/5 
8/7 


0.332 
0.303 
0.342 


02 


55 


5) 


" 




5) 


55 


1 


0.475 


55 


55 


8/3 

8/5 
8/7 


NA'^ 
0.165 
0.203 


03 


55 


55 


55 




55 


55 


1.581 


0.391 


55 


" 


8/3 
8/5 
8/7 


NA 
0.148 
0.171 


04 


55 


55 


55 




5! 


55 


2.5 


0.324 


55 


55 


8/3 
8/5 
8/7 


NA 
0.088 
0.142 


05 


55 


55 


55 




55 


55 


3.952 


0.270 


55 


55 


8/3 
8/5 
8/7 


NA 
0.084 
0.125 


06 


55 


55 


55 




55 


55 


6.248 


0.226 


55 


55 


8/3 
8/5 
8/7 


NA 
NA 
0.093 


07 


55 


55 


55 




5! 


5) 


9.878 


0.190 


55 


55 


8/3 

8/5 
8/7 


NA 
NA 
0.078 


08 


55 


55 


55 




55 


55 


15.62 


0.161 


55 


55 


8/3 
8/5 
8/7 


NA 
NA 
0.062 



* Wavelength of a sampled mode. Other modes exist and were measured that are not listed here. In runs Ola— Olf, the listed mode is the 
strongest mode.'' Runs Ola— Clf experiment with the size and spatial resolution of the simulation box.'^ NA denotes a mode whose Fourier 
amplitude does not grow smoothly and exponentially over the duration of the simulation. In most cases the amplitude either remains 
roughly constant or decreases. 



q = Q disks. Cabot (1984) finds that baroclinic, q = disks arc hncarly unstable to non-axisymmetric perturbations; 
that the "barochnic instabihty draws its energy from radial excursions by fluid elements" (radial excursions made 
possible by the Coriolis force); and that maximum growth rates uii are of order the maximum vertical shearing frequency 
v[\sx\dvy/ dz\ (sec his Tabic 1 and the discussion following his equation 22b). For our dusty layers, Taax.\dVy/dz\ ~ 

fl/^/Ri; support for the inverse square-root dependence on Ri can indeed be found in Figure 4.* 

Knobloch & Spruit (1986) extend the work of Cabot (1984) by restoring the effects of Keplerian shear to baroclinic 
disks. They find that radial shear tends to stabilize the flow to non-axisymmetric disturbances. For their particular 
choice for the form of Vy{z) oc — z^, the linear baroclinic instability is defeated when baroclinity is small {\dvy/dz\ <C O) 
and when the disk is strongly stable to thermal convection in the vertical direction. These stability requirements are 
usually satisfied in disks without gradients in mean molecular weight (such disks are their main concern) . Knobloch & 
Spruit (1986) speculate that instability occurs when baroclinity is large. In the context of their disks, large baroclinity 
demands variations in specific entropy on short radial length scales < h, the vertical height of disk gas. Rapid radial 
variations are also required for instability by Arlt & Urpin (2004), who study axisymmetric perturbations. 

But large baroclinity can also be achieved by having a strong molecular weight gradient in the vertical direction, as 
exists in a highly settled dust layer. We turn now to simulating such highly baroclinic, radially shearing fiows. 

3.2. Coriolis + Radial Shear: q = 3/2 

Having found simulation parameters {Ly,Lz,Ny,Nz) that produce realistic results for g = 0, we use those same 
parameters for our radially shearing g = 3/2 simulations. Since the radial wavemimber of a disturbance grows at 
rate kx = qfloky, our resolution in x should be at least comparable to that in y. Our standard run parameters are 
{Lx, Ly, Lz) = (8, 8, 4)zniax and {Nx, Ny, N^) = (64,64,32). Table 2 lists the various experiments, all of which start 
with Ri = constant flow as described in §2.2. Whether the dust layer turns over within the run duration is indicated 
in the table. 

Our shearing rates are much larger than those considered by Cabot (1984), who analyzed pure gas disks with no gradients in mean 
molecular weight. When = 0, the maximum vertical shear is max\dvy/dz\ = ^Qh/r, where h Cs/Q < r is the hydrostatic disk 
thickness and ^ is a number typically of order 0.1. The smallness of ^ arises because two effects compete and nearly cancel: both the radial 
gravity —d^/dr and the radial pressure acceleration —{l/p)dP/dr decrease in magnitude with increasing height. 
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Fig. 4. — Fourier amplitudes versus time for the \y/zms.yi = 8/7 mode, in simulations of varying Ri. At fixed x = 0.0625zmax, and given 
2 and t, we compute the fast Fourier transform (FFT) of t;^ in the y-direction. What is plotted at a given t is the maximum value of the 
FFT (for the chosen \y) over all z. Growth rates loi are the slopes of lines fitted to the linear growth phase. 
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Fully Shearing [q = 3/2) Simulations 
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Maximum initial perturbation velocity applied to each grid cell. From cell to cell, the perturbation velocity points in a random 
direction, and its magnitude is drawn randomly from a uniform distribution from to the maximum indicated in the table. For 
comparison, the maximum vertical shearing velocity is about iimax ~ 0.05cso-^ For runs in which the layer does not overturn, 
we perform the same Fourier analysis that we do for Coriolis-only runs, and detect no growing modes at all.'^ Increasing /lo 
(equivalently, E^/Sg) at fixed Ri steepens the dust density profile pd{z), whose resolution then requires greater A^'^.'* We do not 
simulate Ri = 0.0625 and Ed/Eg = 0.05, since for such parameters no = 79.5; dust self-gravity, which we do not account for, 
would be significant in this case. 
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Fig. 5. — Snapshots of the (y, z) flow for Ri = 3.952 and /io = 0.270, including only Coriolis forces (q = 0), at times t = (top), t = OSQq ^ 
(middle), and t = 170f2^ (bottom). Greyscale denotes the density of dust only, normalized to the density of dust at the midplane at t = 0. 
Velocity vectors are shown in the frame rotating with dust-free gas (see top panel) . The lone arrow at the top of each panel gives the length 
of "^max for that panel only. Data arc based on a simulation for which (Lx-, Ly., L^) — (0.25, 8. 'l)2^niax 

and {N^,Ny, N^) = (2, 128, 64). 

Results for q = 3/2 differ dramatically from those for g = 0. Whereas instability characterizes all values of Ri for 
g = 0, we find that Ri must fall below a critical value, Ri„it ~ 0.1, for the layer to overturn when g = 3/2 and when 
Ed/Eg = 0.01 (see §3.2.2 for experiments that vary E^/Eg). Figures 6 and 7 demonstrate this point: the flow for 
Ri — 0.0625 eventually mixes whereas that for Ri — 0.125 keeps the dust layer intact for as long as = IQQUq^ . The 
same Fourier analysis of §3.1, applied to the latter run, reveals no growth of any mode. Repeating the Ri = 0.125 run 
at higher resolution — {Nx,Ny,Nz) = (128, 128,64) — confirms these results (run S2a). 



3.2.1. Sensitivity to Initial Velocity Perturbations 

Naturally, reducing the initial velocity perturbations increases the time required for the layer to overturn. The 
dependence is very slight; lowering the maximum perturbation velocity from 10~^Cso to 10~^Cso increases the time to 
overturn by 10-20ri^^ (runs S3a and S4a). (Nonetheless, eliminating the initial perturbations altogether produces no 
evolution in the shearing box whatsoever.) 

The instability grows fastest at ''eo-rotation," where Vy = 0, and spreads radially inward and outward from that 
location. This behavior is evident in Figure 8, and we observe it in all our unstable runs. Upon removing all initial 
velocity perturbations in a narrow annulus surrounding co-rotation, and leaving the perturbations in place everywhere 
else, we find that the flow outside co-rotation still overturns, though more slowly. Eventually the instability spreads 
to co-rotation and mixes dust uniformly throughout the entire box. This experiment assures us that the instability is 
not an artifact of Vy = 0. At present we attribute all this behavior to the fact that our fixed-amplitude initial velocity 
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Fig. 6. — Snapshots of the {y,z) flow, at fixed x = 0, for Ri = 0.0625 and ^to = 1.989, in a fully shearing box {q = 3/2), at times 
t = (top), t = (middle), and t = 37Qq^ (bottom). Greyscale denotes the density of dust only, normalized to the density of 

dust at the midplane at t = 0. Velocity vectors are shown in the frame rotating with dust- free gas (see top panel). The lone arrow at 
the top of each panel gives the length of Vma-x for that panel only. Data are from simulation S3 for which (Lx, Ly, Lz) = (8, 8,4)zmax, 
{Nx,Ny,Nz) = (64,64,32), and Sd/Sg = 0.01. 

perturbations are fractionally larger near co-rotation; i.e., the initial perturbation Svy ^ lO^^Cso is a greater fraction 
of Vy near co-rotation than elsewhere. It would therefore be natural to expect the instability to manifest itself most 
quickly at co-rotation. 



3.2.2. Varying the Height-Integrated Metallicity T,d/T,g 

For solar or super-solar values of ^df^g = 0.01-0.05, the value of Riait changes little from 0.1, as runs S5-S8 of Table 
2 attest. However, iJicrit changes substantially for sub-solar metallicities — it decreases to ~0.02 for Sd/Sg = 0.001 
(runs S9-S12). This result indicates that the Richardson number alone is an inadequate predictor of instability under 
general circumstances. Still, because sub-solar metallicities seem less relevant for planet formation, we have not pursued 
the question of what should replace the Richardson criterion, and we content ourselves with citing Ricrit ~ 0.1 with 
the understanding that this result applies only for solar and moderately super-solar metallicities. 

Values of /Uq and S^/Sg at fixed Ri = 0.1 « Ricin are displayed in Figure 9. In a minimum-mass nebula for which 
Sd/Sg = 0.01, Ri = 0.1 corresponds to a Toomre Q w M/[2nr^pgf){l + /io)] « 25, independent of disk radius if 
PgO oc r~^, as is approximately the case in standard nebular models. Under these conditions, our neglect of the dust 
layer's self-gravity is well justified. Only when Sd/Sg > 0.05 (/zq > 30; see run S8) does dust self-gravity become 
important (Sekiya 1998). 
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Fig. 7. — Snapshots of the {y, z) flow, at fixed x = 0, for Ri 



mox> 

0.125 and fiQ 



1.31, in a fully shearing box {q = 3/2), at times t = (top), 



500(1 (bottom). Greyscale denotes the density of dust only, normalized to the density of dust at the midplane 
at f = 0. Velocity vectors arc shown in the frame rotating with dust-frce gas (see top panel). The lone arrow at the top of each panel gives 



the length of i>max for that panel only. Data are from simulation S2 for which (Lx, Ly,Lz) ■ 



0.01. 



,8,4)2;^ax, {Nx,Ny,N,) = (64,64,32), and 



3.2.3. Baroclinic vs. Kelvin- Helmholtz 



We cannot say whether the observed turnover of the dust layer in a fully shearing disk is better ascribed to the 
Kelvin-Helmholtz instability (as studied by, e.g., Chandrasckhar 1961) or to the baroclinic instability (as studied by 
Cabot 1984, Knobloch & Spruit 1986, and Arlt & Urpin 2004). Both instabilities rely on the vertical shear. For 
Ri w 0.1, the maximum vertical shearing frequency ^Q/y/Ri exceeds, by factors of a few, both the maximum Brunt- 
Vaisala frequency ~r2 (thereby satisfying the traditional criterion for the KHI) and the rotation frequency (thereby 
satisfying the criterion for "large baroclinity," as defined by Knobloch & Spruit 1986). An analytic criterion for the 
baroclinic instability, relevant for non-axisymmetric perturbations, that is valid near the midplane of a radially shearing 
disk and for arbitrary values of the vertical shear, is not known (Knobloch & Spruit 1986). 

4. SUMMARY 

We have performed numerical simulations of flows near the dust-rich midplancs of protoplanetary disks. Dust 
particles are assumed small enough that they are perfectly entrained in gas, and disk self- gravity is neglected. The 
midplane flow is doubly shearing; the rotational velocity varies with height according to the vertical gradient in 
dust density, and it also varies with radius according to the usual Keplerian shear. The simulations are necessarily 
three-dimensional and performed in a shearing box. 

Despite the complications introduced by rotation and tidal gravity, the Richardson criterion for the onset of the 
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Fig. 8. — Snapshots of the midplanc {x,y) flow, at fixed 2 = 0, for Ri = 0.0625 and /.lo = 1.989, in a fully shearing box {q = 3/2), at 
times t = (top), t = 24Qq^ (middle), and t = 37Q^^ (bottom). Grcyscalc denotes the density of dust only, normalized to the density 
of dust at the midplane at t = 0. Velocity vectors arc shown in the frame rotating at Uq. Note how Vy = at x < 0, not at the usual 
X = 0, a, consequence of the imposed background radial pressure gradient [dP/dr)t=o that enforces sub-Keplerian flow. The lone arrow 
at the top of each panel gives the length of Vmax for that panel only. Data are from simulation S3 for which {Lx, Ly, Lz) = (8, 8,4)zmax, 
{Nx,Ny,Nz) = (64,64,32), and S^/Eg =0.01. 

Kelvin-Helmholtz instability still proves useful for solar to moderately super-solar height-integrated metallicities. Dust 
layers characterized by constant Richardson number and S^^/Eg « 0.01-0.05 overturn and mix if Ri < 0.1, but remain 
intact at larger Ri. This result contrasts with the situation when only the Coriolis force is accounted for and the radial 
shear is suppressed: in that case, unstable modes persist for Ri as large as 16. But because these "Coriolis-only" (and 
likely baroclinic; Cabot 1984) modes at large Ri grow at rates that are substantially slower than the Kepler strain rate 
of 3fl/2, they are stabilized and rendered impotent by the Kepler shear. 

In hindsight, our guess (§1) that the critical Ri dividing stability from instability might only change by a factor of 
order unity proved correct, though only for solar to moderately super-solar metallicities: we find that i?icrit for our 
doubly shearing flows is about half that of the traditional value of 1/4, if S^/Sg « 0.01-0.05. This is sensible insofar 
as Riciit ~ 0.1 gives a vertical shearing frequency dvy/dz that is a few times faster than any of the other frequencies 
of the problem: the Brunt- Vaisala frequency for vertical oscillations (this reaches a maximum of ~f2 at 2; ~ Zd), the 
Kepler strain rate (3ri/2), and the Coriolis turning frequency (2^). Deeper insights, including an understanding of 
why Riciit ^ 0.1 for ^d/^g ^ 0.01, might be gained by study of the geophysical literature, where doubly shearing, 
rotating flows are commonplace (see, e.g.. Chapter 7 of Pedlosky 1979, which discusses baroclinic instabilities). 

We conclude that for a standard, height-integrated, solar metallicity of S^/Eg = 0.01, the dust density at the 
midplane can be at most ^l.Sx greater than the gas density. If we assume a gas density appropriate to a minimum- 
mass disk, then the maximum total density of gas and dust is still too low, by a factor of -^25, for the marginally stable 
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Fig. 9. — One-to-one correspondence between the central concentration //o and height-integrated metallicity S^j/Eg at fixed Ri = 0.1 » 
Kicrit- Once the height-integrated metaUicity is several times solar, the midplane dust-to-gas density (1 + /io)PgO can approach Roche 

densities (see §4). 

layer to become gravitationally unstable. The literature discusses two remedies. The first is to find ways of enhancing 
Ed/Eg by factors of 3-10, thereby increasing the dust-to-gas ratio at the midplane by factors of ^30 at fixed Ri = 0.1 
(Sekiya 1998; Youdin & Shu 2002; see also our Figures 3 and 9, and run S8 in Table 2). Possible means of increasing 
the height-integrated metallicity include (i) decreasing through photoevaporation of gas, (ii) increasing E^ by radial 
drift and pile- up of particles (Youdin & Shu 2002; Youdin & Chiang 2004) or (iii) increasing E^ by radiation blow-back 
of grains into the rim of a transitional disk (Chiang & Murray-Clay 2007). Gas giants that form in such a metal- rich 
environment might be expected to have metallicitics enhanced above the solar value by similar factors of 3-10. Indeed 
some hot Jupiters have remarkably metal-rich interiors (Sato et al. 2005). Alternatively, one can relax the assumption 
of perfect coupling between dust and gas, and exploit drag instabilities that result from finite momentum stopping 
times and the backreaction of dust on gas (Youdin & Goodman 2005; Johansen et al. 2007). 

We thank Anders Johansen, Bryan Johnson, Eve Ostriker, Jack Wisdom, and Andrew Youdin for helpful discussions. 
We owe to Jeremy Goodman an insightful referee's report that brought to light the baroclinic instability. This work 
was supported by NSF grant AST-0507805. 
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